home *** CD-ROM | disk | FTP | other *** search
/ NeXT Education Software Sampler 1992 Fall / NeXT Education Software Sampler 1992 Fall.iso / SoundAndMusic / cmix / lib / hplset.c < prev    next >
Encoding:
C/C++ Source or Header  |  1991-12-09  |  4.1 KB  |  168 lines
  1. #include <math.h>
  2. #include "../H/complexf.h"
  3. #include "../H/ugens.h"
  4.  
  5. hplset(xxlp,dur,dynam,plamp,seed,sr,init,q)        
  6. float xxlp,dur,dynam,plamp,seed,sr,*q;
  7. {
  8. //c...............need to replace rand unit    
  9. //c.................compile with order -lF77 -lm    
  10.     float xlp,freq,pi,w,tau,n1,ga1,x,b,s2,s1,rho,c1,pa,pc,cc,c,fl,fu,fm,rl,gl,gl2,w2,r1,r2,r;
  11.            complex j,v1,v2,v3,v4,v5,v6,cexp(),xcmplx(),xadd(),xdivide(),xmultiply(),xsubtract(),smultiply();
  12.     float frand(),zz;
  13.      int n,len,m;
  14.     struct {double x,y;} z; 
  15. //printf("in hplset%f %f %f %f %f %f %d %f\n",xxlp,dur,dynam,plamp,seed,sr,init,q);
  16.            xlp=xxlp;
  17.       freq=1./xlp;   
  18.       pi=3.141592654;
  19.       w=2.*pi*freq;  
  20.       tau=dur/6.91;  
  21.       n1=(int)(sr/freq-.5);
  22. L21:    ga1=exp(-(n1+.5)/(tau*sr));   
  23.       x= 1. - pow(ga1,2.);   
  24.       b=2.*(1.-cos(w/sr));
  25.       s2=b*b-4.*b*x ;   
  26.       if(s2<0) goto L11;
  27.       s1=sqrt(s2)/(2*b);  
  28.       s1=.5+s1 ;
  29.       rho=1.0 ; 
  30.       if(s1<1) goto L10;
  31. L11:   s1=.5;    
  32.       rho=exp(-1./(tau*freq))/ABS(cos(w/(2.*sr)));  
  33. L10:   c=1.-s1;  
  34.       c1=s1;    
  35.       xlp=sr/freq;   
  36.       pa=(-1)*sr/w*atan(-c*sin(w/sr)/((1.-c)+c*cos(w/sr)));  
  37.       pa=1.-pa ;
  38.       n=(int)(xlp-pa);
  39.       if(n == n1) goto L20;
  40.       n1=n;
  41.       goto L21;  
  42. L20:    pc=(sr/freq)-n-pa;  
  43.       cc=sin((w/sr-w/sr*pc)/2.)/sin((w/sr+w/sr*pc)/2.) ;
  44.       fl=20.;   
  45.       fu=sr/2.; 
  46.       fm=sqrt(fl*fu);
  47.       rl=exp(-pi*dynam/sr);                                               
  48.        j=xcmplx(0.,-1.);
  49.        v1 = smultiply(rl,cexp(smultiply(2*pi*fm/sr,j)));
  50.        v2 = xsubtract(xcmplx(1.,0.),v1);
  51.        z.x = v2.re; z.y = v2.im; 
  52.        zz = cabs(z);
  53.        gl = (1.-rl)/zz;
  54.        gl2=gl*gl;
  55.       w2=w/2;  
  56.       r1=(1.-gl2*cos(w/sr))/(1-gl2);
  57.       r2=2*gl*sin(w2/sr)*sqrt(1.-gl2*(pow(cos(w2/sr),2)))/(1.-gl2);
  58.       r=r1+r2; 
  59.       if(r >= 1) r=r1-r2;
  60.       q[1]=n+20;
  61.       len=q[1] ; 
  62.       if(init) {
  63.         for(m=20; m<len; m++) {
  64.                 seed=frand(seed);
  65.             x=seed*2-1.;
  66.                 q[m-1]=plamp;  
  67.                     if(x < 0) q[m-1]= -plamp; 
  68.             }
  69.               q[4]=0.; 
  70.               q[5]=0.; 
  71.               q[6]=0.; 
  72.               q[7]=0.; 
  73.     }
  74.  
  75.       q[1]=q[1]-1;
  76.       if (init) q[0]=q[1];
  77.       q[2]=c ;  
  78.       q[3]=1.-c;
  79.       q[8]=r;  
  80.       q[9]=rho;
  81.       q[10]=cc ;
  82.       q[11]=plamp;   /* for bowplucks */
  83.  
  84.     }
  85. float frand(x)
  86. float x;
  87. {
  88.       int n;
  89.       n=x*1048576.;
  90.       return((float)((1061*n+221589) % 1048576)/1048576.);
  91. }
  92. complex cexp(a)
  93. complex a;
  94. {
  95.     complex b;
  96.     b.re = exp(a.re) * cos(a.im);
  97.     b.im = exp(a.re) * sin(a.im);
  98.     return(b);
  99. }
  100. complex xcmplx(a,b)
  101. float a,b;
  102. {
  103.     complex c;
  104.     cmplx(a,b,c);
  105.     return(c);
  106. }
  107. complex xdivide(a,b)
  108. complex a,b;
  109. {
  110.     complex c,d,e,f;
  111.     conjugate(b,c);
  112.     multiply(a,c,d);
  113.     multiply(b,c,e);
  114.     f.re = d.re/e.re;
  115.     f.im = d.im/e.re;
  116.     /*prtcmplx(a); prtcmplx(b); prtcmplx(c); prtcmplx(d); prtcmplx(e); prtcmplx(f);*/
  117. }
  118. complex xmultiply(a,b)
  119. complex a,b;
  120. {
  121.     complex c;
  122.     c.re = a.re * b.re - a.im * b.im;
  123.     c.im = a.re * b.im + a.im * b.re;
  124.     return(c);
  125. }
  126. complex xadd(x,y)
  127. complex x,y;
  128. {
  129.     complex z;
  130.     z.re = x.re + y.re;
  131.     z.im = x.im + y.im;
  132.     return(z);
  133. }
  134. complex xsubtract(x,y)
  135. complex x,y;
  136. {
  137.     complex z;
  138.     z.re = x.re - y.re;
  139.     z.im = x.im - y.im;
  140.     return(z);
  141. }
  142. complex smultiply(x,y)
  143. float x;
  144. complex y;
  145. {    complex z;
  146.     z.re = y.re * x;
  147.     z.im = y.im * x;
  148.     return(z);
  149. }
  150.  
  151. qnew(freq2,q)                                          
  152. float freq2,*q;
  153. {
  154.     float w,xlp,pa,pc,c;
  155.     int n;
  156.        c=q[2];
  157.       w=2.*PI*freq2 ;                                                    
  158.       xlp=SR/freq2;                                                      
  159.       pa=(-1)*SR/w*atan(-c*sin(w/SR)/((1.-c)+c*cos(w/SR)));              
  160.       n=(int)(xlp-pa);  
  161. printf("in qnew %d %f %f %f %f %f %f %f\n",n,xlp,pa,c,SR,freq2,q[1],q[10]);                                                  
  162.       pa=1.-pa;                                                          
  163.       pc=(SR/freq2)-n-pa;                                                
  164.       q[10]=sin((w/SR-w/SR*pc)/2.)/sin((w/SR+w/SR*pc)/2.);                 
  165.       q[1]=n+19;            
  166. printf("in qnew after changes %f %f\n",q[1],q[10]);                                             
  167. }                                                           
  168.